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Abstract 

Composite  laminates  could  be  made  to  diffuse  penetration  energy  at  predetermined  rates  in  addition  to  locations. 
Materials  with  different  hardness  and  toughness  can  be  stacked  in  the  appropriate  sequence  to  create  a  more 
widespread  damage  pattern.  The  kinetic  energy  per  unit  volume  of  the  laminate  could  then  be  reduced  to  prevent 
complete  perforation. 

Analyzed  in  detail  is  the  failure  sequence  of  a  tungsten  rod  (47  mm  long  and  14  mm  in  diameter)  impacting  a 
three-layered  composite  laminate  at  1100  m/s.  The  top  layer  is  a  6  mm  thick  ceramic  (A1203).  Sandwiched  in 
between  is  a  10  mm  thick  rubber  layer  backed  up  by  a  4  mm  thick  4340  steel  plate.  At  approximately  10  ps,  the 
tungsten  rod  is  already  badly  shattered,  whereas  the  ceramic  is  only  fragmented  in  a  layer  next  to  the  surface.  The 
ceramic  reflects  the  impact  waves  which  cause  damage  in  the  tungsten  rod.  Within  the  next  1.50  ps,  i.e.  21.5  ps,  the 
dissipated  energy  density  in  the  rubber  increased  two  orders  of  magnitude  from  32.11  to  33.25  X  10 2  Pa  while  those 
in  the  ceramic  and  steel  are  much  smaller.  The  highest  temperature  in  the  ceramic  reached  about  555°C  and  that  in 
the  back-up  steel  plate  is  only  275°C.  Upon  complete  fragmentation  of  the  tungsten  rod,  a  through  the  thickness 
cavity  is  predicted  in  the  ceramic  while  a  mushroom  head  shaped  crack  is  predicted  in  the  rubber.  Interfacial  cracks 
are  widespread  but  the  steel  layer  remained  intact. 

Presentation  is  divided  into  two  parts.  Part  I  describes  the  basic  scheme  of  the  methodology  and  the  nonhomoge- 
neous  character  of  response  within  each  layer  of  the  composite  laminate.  That  is  the  constitutive  relation  for  each 
material  element  is  derived  according  to  the  local  strain  and  strain  rate  rather  than  preassumed.  Thermal/mechani¬ 
cal  interactions  are  synchronized  by  application  of  the  isoenergy  density  theory  in  nonequilibrium  thermomechanics. 
At  0.15  ps,  concentration  of  energy  is  detected  at  the  impact  site  and  ceramic/rubber  interface.  The  finite  element 
grid  pattern  had  to  be  remeshed  after  50  ps  to  account  for  the  discontinuities  created  by  damage  at  the  interphases. 
Modelled  in  Part  II  are  the  evolution  of  damage  in  steps  until  the  tungsten  rod  is  completely  destroyed. 
Nonequilibrium  data  such  as  load  strain  rates,  temperatures,  dissipation  energies,  etc.  are  given  numerically  and 
displayed  graphically. 


1.  Introduction 

Dynamic  response  of  composite  laminates  cov¬ 
ers  an  overwhelmingly  wide  range  of  problems 
that  could  vary  from  the  rudimentary  tests  to  the 


sophisticated  analyses.  Many  of  the  past  works 
can  be  found  in  [1-3].  Assessment  of  the  nonho- 
mogeneous  character  of  the  damage  process  be¬ 
comes  problematic  because  the  material  constitu¬ 
tive  relation  and  failure  criterion  must  be  as- 
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sumed  independently  in  the  classical  continuum 
mechanics  approach.  There  is  no  rational  means 
to  account  for  the  disparities.  In  principle,  pre¬ 
dictive  models  should  have  the  capability  to  scale 
the  base  reference  material  data  to  situations 
other  than  those  tested.  That  is  the  conversion  of 
low  strain  rate  uniaxial  data  to  high  strain  rates 
and  low  temperature  data  to  high  temperature 
levels.  Such  a  scheme  is  inherent  in  the  isoenergy 
density  theory  [4,5]. 

Damage  patterns  for  a  single  4340  steel  plate 
20  mm  thick  impacted  by  a  tungsten  rod  at  a 
speed  of  1100  m/s  have  been  obtained  [6].  Using 
only  tensile  specimen  data  at  ambient  condition, 
local  temperature  of  640°C  and  latent  heat  (strain 
rate  of  energy  dissipation  density)  of  11  MPa 
were  predicted  at  0.35  p,s.  Perforation  of  the 
plate  by  plugging  was  predicted  in  addition  to  the 
phase  transformation  of  the  material  in  a  narrow 
band  around  the  plug.  Thickness  of  the  plug  is 
about  40%  of  that  for  the  plate  while  a  33% 
reduction  of  the  tungsten  rod  was  estimated.  The 
rod  velocity  decreased  to  760  m/s  at  0.30  |xs. 
Such  details  have  also  been  made  available  for 
the  hypervelocity  impact  of  aluminum  plate  [7]. 
Local  melting  and  vaporization  of  the  aluminum 
were  predicted  where  the  estimated  strain  rates 
reached  as  high  as  108-109  s-’.  These  findings 
are  in  general  agreement  with  the  experimental 
findings.  Validation  of  nonequilibrium  tempera¬ 
tures  with  experiments  can  be  found  in  [5]  for 
uniaxial  specimens  under  static  loadings.  More 
recently,  dissipation  energy  densities  were  mea¬ 
sured  [8]  for  characterizing  the  failure  behavior  of 
composite  materials.  They  describe  the  degree  of 
internal  damage  in  more  realistic  terms. 

Presented  is  Part  I  of  this  study  concerning  the 
impact  damage  of  a  tungsten  rod  striking  a 
three-layered  ceramic/rubber/steel  composite 
laminate.  The  stress  and  strain  response  for  each 
material  element  is  determined  from  a  knowledge 
of  the  rate  change  of  local  volume  to  surface 
which  couples  the  exchange  between  surface  and 
volume  energy.  It  also  synchronizes  the 
thermal/mechanical  disturbances  caused  by  the 
impact.  The  IDA  (isoenergy  density  analysis)  pro¬ 
gram  [9]  is  used  to  formulate  the  problem  for 
analyzing  the  nonhomogenous  response  of  the 


composite  laminate.  The  sequence  of  the  damage 
process  will  be  discussed  in  Part  II  [10]  where  the 
calculations  are  carried  out  until  the  analysis 
predicts  complete  destruction  of  the  tungsten  rod. 


2.  Isoenergy  energy  theory:  finite  element  size 

Two  distinct  features  of  the  isoenergy  density 
theory  are  [4,5]: 

•  Element  size  remains  finite  and  is  determined 
from  the  condition  of  isoenergy  density  state. 

•  Temperature  is  obtained  directly  from  the  rate 
change  of  volume  with  surface,  isostrain  and 
dissipation  energy  density. 

The  term  isoenergy  refers  to  the  situation  where 
the  same  energy  is  transmitted  across  the  orthog¬ 
onal  surfaces  of  a  continuum  element.  Energy 
density  could  still  vary  from  element  to  element. 
What  the  isoenergy  density  theory  establishes  is 
the  unique  correspondence  of  energy  states  re¬ 
gardless  of  their  stress  and/or  strain  dimension¬ 
ality.  No  other  continuum  mechanics  theories  with 
energy  dissipation  possess  such  capability.  The 
classical  theory  of  plasticity  leaves  out  the  dilata- 
tional  portion  of  the  energy  density  when  making 
correlation  between  the  uniaxial  and  multiaxial 
stress-strain  states.  Separation  of  dilatation  and 
distortion  is  borne  out  from  the  property  of  linear 
elasticity;  it  is  not  valid  for  nonlinear  theories. 

2.1.  Nonlocal  effect 

In  the  isoenergy  density  theory,  the  traction  7] 
on  a  tetrahedron  referred  to  the  rectangular  co¬ 
ordinates  (x,  y,  z)  is  given  by 

fi.rllnl+P(u,-hl)^'):  (1) 

where  a-  are  stress  components,  ii,  are  the  accel¬ 
erations  and  ht  the  body  force  components  with 
p  being  the  mass  density.  The  displacement  com¬ 
ponents  are  n;.  Eq.  (1)  differs  from  classical  con¬ 
tinuum  mechanics  in  that  the  rate  changes  of 
volume  with  surface  (AE/A,4)„  or  2^  are  no 
longer  assumed  to  vanish  in  the  limit.  It  provides 
the  additional  length  scale  parameter  that  ac- 
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counts  for  nonlocal  effects.  The  interaction  of  the 
inertia  term  pii,  with  T~n  becomes  crucial  at  the 
impact  site  where  surface  energy  is  being  con¬ 
verted  into  volume  energy. 

2.2.  Energy  transmission  across  surface 

When  two  bodies  are  impacted,  energy  trans¬ 
fers  from  one  body  to  the  other  across  the  con¬ 
tact  surface  or  interface.  There  is  an  exchange  of 
surface  and  volume  energy  density.  Let  the  for¬ 
mer  be  denoted  by  ( AW/AA ),  or  ^  and  the 
latter  by  AiV/AV  or  W.  They  are  related  as 

=  i  =x,  y,  z.  (2) 

The  quantities  represent  the  energy  prevail¬ 
ing  on  the  three  orthogonal  surfaces  of  an  ele¬ 
ment,  while  2^  are  the  corresponding  compo¬ 
nents  of  the  vector  for  the  rate  change  of  volume 
with  surface.  For  an  isoenergy  element  the  three 
components  ^  (i  =  1,  2,  3)  are  required  to  be 
the  same.  The  significance  of  this  requirement 
will  be  elaborated  subsequently. 

2.3.  Isoenergy  elements 

The  size  and  orientation  of  an  isoenergy  ele¬ 
ment  can  be  determined  by  letting 
&i=*S’2*=S’3-*S,,0,  i  =  1,  2,  3,  (3) 

that  is  the  same  y0  would  prevail  on  all  of  the 
three  orthogonal  surfaces.  Since  W  is  a  scalar,  it 
follows  from  Eq.  (2)  that 

The  rate  change  of  volume  with  surface  for  an 
isoenergy  element  is  also  the  same  in  all  direc¬ 
tions.  These  different  size  elements  could  be  con¬ 
nected  to  form  a  continuous  body  such  that  upon 
deformation  there  will  be  no  gaps  between  the 
conceptual  interfaces  of  the  adjacent  elements. 
The  proof  for  the  existence  of  the  isoenergy  sur¬ 
face  made  of  these  elements  is  given  in  [4]. 

2.4.  Field  equations 

In  order  to  emphasize  the  nonlocal  character 
of  the  isoenergy  density  theory,  Eq.  (1)  can  be 


rewritten  in  terms  of  the  isostress  components  t,  • 
as 

f,  =  Tj,nj,  Tij  #  Tji,  (5) 

in  which 

Tp  ~  +  p(w,  —  hf)^,  (6) 

and  r  is  not  symmetric.  The  equations  of  equilib¬ 
rium  expressed  in  terms  of  Ty  have  the  same 
formal  appearance  as  those  in  the  classical  theo¬ 
ries,  i.e. 

3  Tji 

—  +phi  =  pui,  (7) 

Hj 

where  £,■  ( i  =  1,  2,  3)  refer  to  the  current  state 
coordinates  of  an  element  as  distinguished  from 
the  reference  coordinates  xl  (i  =  1,  2,  3).  The 
conservation  of  energy  requires  that 

j  phjtij  dV  +  j  TjUf  d  A 

,  d  , 

=  J  piiiuidV+—J  W  dV,  (8) 

in  which  A  is  the  volume  bounded  by  the  surface 
I  and 


Note  that  the  volume  energy  density  function  W 
is  not  an  elastic  potential.  It  applies  to  an  irre¬ 
versible  and  dissipative  process.  The  dot  repre¬ 
sents  differentiation  with  time.  The  nine  displace¬ 
ment  gradients  are  etj  given  by 

3  U: 

eU  =  TT ’  L  j=  1,  2,  3.  (10) 

No  limitations  have  been  placed  on  the  finiteness 
and/or  magnitude  of  the  deformation.  When  W 
is  interpreted  as  an  isoenergy  density  function,  it 
has  the  unique  property  that 

W=  Jtu  den  =  Jr, 2  de12  = 

=  /t23  de23^>  fr  de.  (11) 

This  means  that  any  one  of  the  nine  pairs 
(rn,  en),  (r12,  en),  ...,(t23,  e23)  or  (t,  e)  can  be 
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used  to  yield  the  same  W.  Hence,  the  uniaxial 
isostress  r  versus  the  isostrain  e  plot  provides  a 
general  description  of  the  energy  state,  even 
though  the  other  stress  components  r,7  are  also 
acting  on  the  element. 

2.5.  Synchronization  of  temperature  and  deforma- 


The  nonequilibrium  temperature  0  can  be 
determined  directly  from  e  and  S:  without 
introducing  the  concept  of  heat.  Here,  3>  stands 
for  the  dissipation  energy  density 
d.&  =  d W-  d.j /,  &  >  0.  ( 12) 

with  stf  being  the  available  energy  density.  The 
dissipated  and  available  energy  are  mutually  ex¬ 
clusive  and  hence  there  is  no  loss  in  generality  in 
the  linear  dependence  state  in  Eq.  (12).  Derived 
in  [4]  is  the  relation 
A0  Use 

<l3> 

where  0  represents  the  nonequilibrium  tempera¬ 
ture  in  contrast  to  T  in  classical  thermodynamics 
that  applies  only  to  equilibrium  states.  The  pa¬ 
rameter  A  relates  the  slope  of  the  isostress  r 
isostrain  e  curve  to  the  rate  change  of  volume  to 
surface  dV/dA  or  7/\  i.e., 

r  =  f\7'de.  (14) 

2.6.  Nonhomogenous  response 

Since  both  A  and  "V  can  be  computed  from  a 
knowledge  of  the  deformation  of  each  isoenergy 
element,  only  the  values  of  r  and  e  or  2^  and  e 
for  the  initial  or  reference  state  need  to  be  known. 
Once  the  loading  steps  are  specified,  adjustment 
on  (t,  e)  can  be  made  for  each  of  the  subsequent 
steps.  In  other  words,  the  reference  state  from  0 
to  1  in  Fig.  1  or  (^°,  e°)  is  assigned  while  the 
subsequent  steps  (2^',  e1),  (2^2,  e2\  etc.,  are  de¬ 
termined,  the  loci  of  which  gives  the  response  of 
a  typical  element.  It  is  a  self-adaptive  procedure 
that  accounts  for  the  inhomogeneous  response  of 
the  system. 


Fig.  I.  Path  of  typical  element  in  isostress  and  isostrain 
domain. 


The  foregoing  procedure  has  been  used  to 
generate  high  strain  rates  uniaxial  data  beyond 
the  range  of  present-day  test  machines.  Data  on 
complete  stress  and  strain  histories  can  be  moni¬ 
tored  only  for  strain  rates  no  higher  than  10  s_1. 
Local  strain  rates  for  the  impact  problem  at  hand 
reach  as  high  as  104  s_l.  Experimental  checks 
can  only  be  made  at  discrete  points  of  the  stress 
and  strain  curve. 


3.  Method  of  approach 

Shown  in  Fig.  2  is  a  three-layer  composite 
laminate  with  one-half  symmetry.  It  is  20  mm  in 
total  thickness  consisting  of  three  different  mate¬ 
rials.  A  ceramic  known  as  corundum  (A1203) 
being  the  top  layer  is  6  mm  thick.  The  middle 
rubber  layer  is  10  mm  thick  while  the  third  layer 
is  made  of  4340  steel  4  mm  thick.  The  staking 
sequence  is  selected  so  that  the  rubber  can  ab¬ 
sorb  and  diffuse  the  initial  surge  in  impact  energy 
while  the  ceramic  being  hard  momentarily  slows 
down  and  breaks  up  an  impacting  tungsten  rod 
with  an  initial  mass  m0  of  0.128  kg  and  velocity 
1 0  of  1100  m/s. 
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Fig.  2.  Composite  laminate  subjected  to  normal  impact. 


3.1.  Material  data  bank 

As  mentioned  earlier,  only  the  initial  mechani¬ 
cal  and  fracture  properties  of  the  composite  lami¬ 
nate  need  to  be  known.  They  are  shown  in  Table 
1.  Once  the  system  is  disturbed  by  impact,  the 
response  of  each  material  element  will  be  derived 
individually  in  accordance  with  the  displacement 
gradients  and  rate  change  of  volume  with  surface. 
The  data  in  Table  1  are  referred  to  as  the  base 
properties. 


Fig.  3.  Material  data  bank  for  tungsten. 


Starting  with  the  base  material  property  of  the 
tungsten  rod,  a  material  data  bank  may  be  con¬ 
structed  such  that  all  points  in  the  (r,  e)  plane 
can  be  located  in  terms  of  'V.  Some  typical 
curves  are  displayed  in  Fig.  3.  A  great  deal  of 
accuracy  can  be  achieved  in  determining  a  given 
state  (V,  e)  for  each  time  increment  by  applica¬ 
tion  of  the  floating  data  bank  scheme.  The  termi¬ 
nal  points  of  all  the  curves  are  governed  by  the 
yield  strength  ays  and  fracture  toughness  Wc  re¬ 
lation  shown  in  Fig.  4  for  the  tungsten.  The  path 
of  a  typical  element  in  the  isoenergy  plane  is 
shown  in  Fig.  3.  It  need  not  follow  any  one  of  the 
curves  in  the  data  bank  which  serves  as  a  visual 
aid  and  to  identify  (r,  e )  with  (‘V,  e). 


Table  1 

Mechanical  and  fracture  base  properties  for  the  three-layer  composite  laminate 


Material  type 

Young’s  modulus 

E  X  1011  (Pa) 

Mass  density 
p(kg/m3) 

Yield  strength 
<rys  X  108  (Pa) 

Ultimate  strength 
<ru  X  108  (Pa) 

Volume  energy  density 
iTc  X  105  (Pa) 

Tungsten 

3.450 

17600 

0.1062 

0.1347 

1.503  x  10~2 

Ceramic 

4.650 

19350 

0.3535 

0.3535 

1.344  X  10-2 

Rubber 

4.0  X  10~7 

1100 

7.500  X  10“6 

1.500  X  10  3 

2.176  X  10“ 1 

4340  steel 

2.068 

7850 

9.70 

11.00 

115.60 
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Critical  volume  energy  density^  (MPa I 
Fig.  4.  Relation  between  yield  strength  and  fracture  toughness 
for  tungsten. 


The  material  bank  for  the  ceramic  layer  is 
shown  in  Fig.  5.  It  consists  of  a  group  of  straight 
lines  because  it  is  very  hard  in  comparison  with 
the  other  materials.  This  does  not  imply  that  all 
elements  in  the  ceramic  layer  follow  a  linear 
stress  and  strain  response.  Linearity  may  not  be 
preserved,  as  the  actual  behavior  of  a  given  ele¬ 
ment  can  cut  across  the  family  of  straight  lines  in 
a  nonlinear  fashion  as  illustrated  in  Fig.  5.  Again, 
the  terminal  points  of  the  stress  and  strain  curves 
in  Fig.  5  are  determined  from  the  Wc  versus  <rys 


Yield  strength  x  106  (Pa) 


Fig.  6.  Relationship  between  fracture  toughness  and  yield 
strength  for  ceramic  in  material  data  bank. 


relation  in  Fig.  6  for  the  ceramic.  The  response  of 
the  rubber  to  load  is  very  different  from  that  of 
the  tungsten  and  ceramic  materials  discussed  pre¬ 
viously.  Fig.  7  shows  that  T’  or  dr/de  can  in¬ 
crease  very  rapidly.  This  feature  is  desirable  from 
the  energy  absorption  rate  viewpoint,  which  will 
become  evident  when  numerical  results  are  made 


Fig.  5.  Material  data  bank  for  the  ceramic  layer. 


Fig.  7.  Material  data  bank  for  rubber  layer. 
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Yield  strength  x  10!  (Pa) 


Fig.  8.  Relationship  between  fracture  toughness  and  yield 
strength  for  rubber  in  data  bank. 


available  subsequently.  A  plot  of  against  o-ys 
for  the  rubber  is  given  in  Fig.  8. 

The  back-up  layer  is  made  of  4340  steel  that 
contains  the  rubber  and  absorbs  the  excess  en¬ 
ergy  transmitted  through  the  first  two  layers.  Fig. 
9  and  10  provide  the  necessary  data  bank  for  the 
expediency  of  the  numerical  calculations. 


Fig.  9.  Material  data  bank  for  4340  steel. 


3.2.  Axisymmetric  deformation 

The  tungsten  rod  is  47  mm  in  length  and  has  a 
diameter  of  14  mm.  Axisymmetry  is  assumed  to 
prevail  at  normal  impact  such  that  the  configura¬ 
tion  in  Fig.  2  should  be  viewed  as  a  body  of 
revolution.  Because  of  axisymmetricity,  only  one 
angle  j8  is  require  to  locate  the  position  of  the 
isoenergy  density  element  in  the  (£,  17,  £)  plane 
as  distinguished  from  the  physical  plane  ( x ,  y,  z) 
or  (r,  0,  z ).  According  to  [4],  the  corresponding 
displacements  are 

«f  =  Wf(£,  uv  =  uv(r1,  C),  uf  =  uf(£,  0, 
(15) 

while  the  five  nonvanishing  deformation  gradi¬ 
ents  are  given  by 

3  Ug  duv  du(  3  Uf 

ef=_3^’  6'n  =  ~d^'  ~  ’  ^=~3 J' 


Critical  volume  energy  density^  (MPa) 


Fig.  10.  Relation  between  yield  strength  and  fracture  tough¬ 
ness  for  4340  steel. 
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The  expressions  for  (dV/dA)l  or  7]  take  the 
form 

r,  =  {2ev +  2(ef  +  e£)(l+ev)  cos  fi 

+  2e,ec  cos2/3  -  (g2  +f[)  sin2/3  -  2 f(g() 
x|[2eI)(l  +e(  cos  p)  +  2e(  cos  p 
-//  sin2/S]  cos  p}  , 

=  {2e„  +  2(ef +  ef)(l  +  e„)  cos  j8 

+  2e(e(  cos2/3  -  (g|  +//)  sin2/3  -  2/fgt} 
x/2[ef  +  ef(l  +e^  cos  j8)]  cos  p 
~  (s|  +//)  sin2)?  -  2/fg J  ,  ( 17) 

=  {2ev  +  2(ef  +  <?f  )(1  +  e^)  cos  p 

+  2e(ec  cos2j8  -  (g|  +/f2)  sin2/3  -  2/fgJ 
X  | [2eTJ(  1  +  e(  cos  P)  +  2ef  cos  0 
-gf  sin2/S  ]  cos  p}  . 

Even  though  the  numerators  and  denominators 
of  Eq.  (17)  contain  only  second-degree  terms  in 
e {,  ev,  etc.,  displacement  gradient  terms  of  order 
higher  than  second  are  included  in  the  expan¬ 
sions  for  5X(i  =  £,  r],  £).  Because  of  axisymmetry, 
the  condition 

2X  =  k  =  const .  (18) 

may  be  assumed  which,  when  applied  to  the 
second  of  Eqs.  (17),  renders 

ev  =  {(K-l)[2(ef+e()  cos  p 

+  2efec  cos2/3  -  2 f(gf  -  (//  +  gf2)  sin2#]} 

x{2[1  +  (ef  +  e()  cos  i8]}  '•  (19) 

Hence,  and  e(  are  independent  since  e v  can 
be  eliminated  in  view  of  Eq.  (19).  The  condition 
of  isoenergy  in  Eq.  (4)  may  be  further  invoked  to 
give 

2(l+ev)(es-ec)  cos  #-(g;-/f2)  sin2#  =  0. 

(20) 


Eq.  (20)  solves  for  j8  or  the  position  of  the  isoen¬ 
ergy  density  element,  since  e e(,  /.-  and  gf  are 
known  functions  of  #,  i.e., 
e(  =  er  cos2 #  +  e:  sin2/3  +  (/,  +  gr)  sin  #  cos  #, 
e(  =  er  sin2#  +  ez  cos2/8  -  ( /.  +  gr )  sin  #  cos  #, 

/.  =  (e:  -  er)  sin  #  cos  #  +/.  cos 2#  -gr  sin2#, 
(21) 

gi  =  (e,~  er )  sin  #  cos  ft  -/,  sin2#  +  gr  cos2#. 
The  quantities  cr,  e,,  /,  and  gr  are  computed 
numerically  by  a  finite  element  procedure  [9]. 

3.3.  Computational  procedure 

The  IDA  program  [9]  is  a  finite  element  for¬ 
mulation  of  the  isoenergy  density  theory.  It  solves 
for  the  displacements  that  are  up-dated  for 
each  increment  of  loading.  The  following  steps 
are  to  be  executed: 

•  Without  loss  in  generality,  start  with  an  initial 
value  of  =  1  mm  and  dr/de  from  the  base 
properties  of  the  material.  This  also  deter¬ 
mines  A. 

•  The  displacements  m,  are  found  to  yield  er,  e,, 
f.  and  gr  in  Eqs.  (21). 

•  Since  es,  ex,  /,  and  g^  can  now  be  related  to 
er,  ex,  fz  and  gr,  the  angle  #  in  Eq.  (20)  can 
be  determined. 

•  This  gives  dV/d  A  or  TC  together  with  e  to 
establish  a  new  (r,  e)  state  by  using  the  initial 
A. 

•  A  segment  of  the  path  from  0  *  to  1  as  indi¬ 
cated  in  Fig.  1  is  established. 

•  The  foregoing  process  can  thus  be  repeated  to 
yield  the  subsequent  segments. 

The  displacement  gradients  er,  ez,  /2  and  gr 
referred  to  the  physical  coordinates  (r,  z)  can  be 
obtained  from  those  referred  to  (£,  £)  using  the 
relations: 

er  =  ef  cos2#  +  ec  sin2#  —  (f(  +gf)  sin2#  cos  #, 
e:  =  e(  sin2#  +  ee  cos2#  +  (/f  +  gf)  sin  p  cos  p, 

f.  =  -  (ef  -  ef)  sin  p  cos  p  +  f(  cos2#  -  gf  sin2#, 

(22) 

g, =  -(e{-ef)  sin  p  cos  P~f(  sin2#  +g^  cos2#. 
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The  displacement  gradients  er  and  ez  should  not 
be  taken  as  the  strain  components  in  the  classical 
continuum  mechanics  theory. 

3.4.  Damage  thresholds 

From  a  knowledge  of  the  data  bank  estab¬ 
lished  earlier,  critical  values  of  y  and  W  or  yc 
and  arc  can  be  computed  for  each  element;  they 
would  change  with  time.  Material  damage  states 
can  be  completely  characterized  by  the  trade-off 
between  'V  and  W. 

Dominant  surface  energy  density.  When  the 
surface  energy  density  exceeds  critical,  i.e.,  y> 
yc,  the  following  two  conditions  prevail: 

•  W-  —  This  corresponds  to  failure  by  frac¬ 


ture  over  a  region,  the  inside  of  which  is  also 
fragmented.  Sliding  nodes  would  be  intro¬ 
duced  between  the  fragments. 

•  W<  Wc  —  Gaps  are  created  by  fracture  be¬ 
tween  material  elements  and  accounted  for  by 
updating  the  grid  pattern. 

Dominant  volume  energy  density.  In  situations 
where  W>WC  the  surface  energy  density  should 
be  compared  with  its  critical  value  yc: 

•  y=*yc  —  Failure  by  massive  fragmentation 
would  be  modelled.  Material  elements  are  dis¬ 
lodged  and  mass  redistribution  to  neighboring 
nodes  are  made. 

•  y  <yc  —  Fragmentation  is  not  as  severe,  but 
mass  redistribution  would  be  needed  to  ac¬ 
count  for  local  damage. 
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Fig.  12.  Element  numbers  for  ceramic /rubber  and  rubber /steel  interface. 


The  initial  time  interval  is  selected  to  be  suffi¬ 
ciently  small  such  that  none  of  the  foregoing 
damage  thresholds  are  exceeded. 


4.  Discussion  of  results:  initial  impact 

The  selection  of  the  initial  time  interval  in¬ 
volves  several  trial  runs  in  order  to  gain  a  knowl¬ 
edge  of  response  sensitivity  for  the  disturbance  in 
the  composite  laminate  system.  Unlike  the  single 
layer  armor  system  [6],  the  changes  in  damage 
pattern  are  more  tolerant  to  time  variations.  This 
is  mainly  due  to  the  damping  effect  provided  by 
the  rubber  layer.  Only  the  results  for  t  =  0.150  ps 
will  be  discussed  in  Part  I  of  the  work.  Presented 


are  the  surface  and  volume  energy  density  con¬ 
tours.  Remeshing  of  the  finite  element  grid  pat¬ 
tern  becomes  necessary  as  the  damage  progresses 
and  alters  the  connectivity  of  the  elements.  The 
latter  will  be  considered  in  Part  II  [10]. 

4. 1.  Contact  at  impact 

The  twelve  node  isoparametric  finite  elements 
embedded  with  Gaussian  points  are  employed  to 
discretize  the  composite  laminate  system,  as  illus¬ 
trated  in  Fig.  11.  The  center  of  axisymmetry  coin¬ 
cides  with  the  sides  of  elements  39,  36,  33,  32  and 
1  such  that  the  adjacent  elements  54,  53,  etc.,  are 
mapped  across  the  center  line  to  retain  symmetry 
in  the  solution.  Special  attention  is  given  to  the 


Table  2 

Isostress,  isostrain  and  isoenergy  density  at  t  =  0.15  p.s 


Material  Element  Isostress  Isostrain 

type  number  r  (MPa)  e  x  l(r5  (m/m) 


Ceramic  23 

24 

Rubber  20 

27 

32 

Steel  3 


-0.00934  -0.00100 

0.00214  0.00024 

0.00695  0.00075 

-0.01372  -0.2368 

-0.01473  -  0.2532 

-0.01091  -0.1923 

-0.2399  -  0.05805 

-0.2196  -0.05301 

-0.2256  -  0.05459 


Surface  energy  density 
■V  (N/m) 

0.00513 

0.00268 

0.00084 

0.53433 

0.61874 

0.29832 

0.03936 

0.04014 

0.03751 


Volume  energy  density 
ST  (Pa) 


0.00116 

0.00060 

0.00018 

0.1943 

0.2241 

0.1086 

0.01984 

0.02004 

0.01869 
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interface  elements  in  Fig.  12  which  are  labelled 
51/48,  50/49,  34/33,  25/26,  22/21,  13/14  and 
10/9  for  the  ceramic/rubber  interface  and  44/43, 
45/42,  31/30,  28/29,  19/18,  16/17  and  7/6  for 
the  rubber /steel  interface.  The  surfaces  between 
the  dissimilar  materials  can  move  relative  to  one 
another. 

Contact  at  impact  is  assumed  to  take  place  at 
0.15  ps.  Results  obtained  from  the  IDA  program 
show  that  the  disturbance  would  first  reach  the 
elements  under  the  impactor.  They  correspond  to 
elements  1,  2,  3  and  4  for  the  steel;  15,  20,  27  and 
32  for  the  rubber;  and  17,  23,  24  and  35  for  the 
ceramic.  Their  relative  positions  are  shown  in 
Fig.  11.  Summarized  in  Table  2  are  the  isostress 
t,  isostrain  e,  surface  isoenergy  density  5?  and 


volume  isoenergy  density  W.  Even  though  their 
amplitudes  are  small,  it  is  worthwhile  to  observe 
two  important  features  of  the  solution  at  initial 
impact.  While  the  rubber  and  steel  are  still  in 
compression,  reflective  tensile  waves  have  already 
reached  elements  24  and  35  in  the  ceramic.  These 
elements  are  directly  under  the  impacting  tung¬ 
sten  rod  whereas  element  23,  being  slightly  away 
to  the  side  of  contact,  is  not  affected  by  wave 
reflection.  The  larger  tensile  stress  is  in  element 
35,  which  is  closest  to  the  interior  or  farthest 
away  form  the  edge  of  the  impactor.  Of  interest  is 
that  the  surface  and  volume  energy  density  in  the 
rubber  are  two  orders  and  one  order  of  magni¬ 
tude  larger  than  those  in  the  ceramic  and  steel, 
respectively.  This  is  a  desirable  feature  of  the 


Fig.  13.  Surface  energy  density  in  laminate. 
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Table  3  Table  4 

Surface  isoenergy  density  values  for  contours  in  Fig.  13  Volume  isoenergy  density  values  for  contours  in  Fig.  14 


Contour 

No. 

S* 

(N/m) 

Contour 

No. 

(N/m) 

Contour 

No. 

sr 

(Pa) 

Contour 

No. 

V 

(Pa) 

1 

0.56X  10~" 

11 

0.75X10  1 

1 

0.18X  10~7 

11 

0.87X10 

2 

0.75X  10~2 

12 

0.82 

2 

0.87X10" 

12 

0.96 

3 

0.15  X 10" 1 

13 

0.90 

0.17x10 

13 

0.10X102 

4 

0.22 

14 

0.97 

4 

0.26 

14 

0.11 

5 

0.30 

15 

0.10X10" 

5 

0.35 

15 

0.12 

6 

0.37 

16 

0.11 

6 

0.44 

16 

0.13 

7 

0.45 

17 

0.12 

7 

0.52 

17 

0.14 

8 

0.52 

18 

0.13 

8 

0.61 

18 

0.15 

9 

0.60 

19 

0.13 

y 

0.70 

19 

0.16 

10 

0.67 

20 

0.14 

10 

0.79 

20 

0.17 

21 

0.15 

21 

0.17 

Fig.  14.  Volume  energy  density  in  laminate. 
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Fig.  15.  Surface  energy  density  near  tungsten-ceramic  inter¬ 


laminate  composite  for  otherwise  the  energy  in 
the  top  and  bottom  layer  would  have  been  higher 
had  the  plate  been  made  of  a  single  material. 

4.2.  Isoenergy  density  contours 

Plotted  in  Figs.  13  and  14  are,  respectively,  the 
S0  and  W  contours  in  the  laminate  composite 
0.15  p.s  after  impact.  Note  that  energy  has  al¬ 
ready  been  transmitted  to  all  parts  of  the  lami¬ 
nate,  although  the  values  of  5?  and  W  are  rela¬ 
tively  low  as  shown  in  Tables  3  and  4. 

Surface  and  volume  energy  concentrations  at 
the  reentrant  corner  where  the  tungsten  rod 
comes  into  contact  with  the  ceramic  layer  are 
illustrated  in  Figs.  15  and  16.  The  2?  and  W 
contours  are  very  closely  packed  at  the  corner. 
Table  5  show  a  difference  of  five  orders  of  magni- 


Fig.  16.  Volume  energy  density  near  tungsten-ceramic  inter¬ 


lude  for  S  at  the  corner  in  comparison  with  the 
surface  energy  density  values  a  distance  away. 

A  similar  situation  prevails  for  the  volume 


Table  5 

Surface  isoenergy  density  values  for  contours  in  Fig.  15 


Contour 

No. 

& 

(N/m) 

Contour 

No. 

se 

(N/m) 

1 

0.19X10  1 

5  12 

0.36 

2 

0.32X10- 

3  13 

0.39 

3 

0.65 

14 

0.42 

4 

0.97 

15 

0.45 

5 

0.13X10 

2  16 

0.49 

6 

0.16 

17 

0.52 

7 

0.19 

18 

0.55 

8 

0.23 

19 

0.58 

9 

0.26 

20 

0.62 

10 

0.29 

21 

0.65 

11 

0.32 
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energy  density.  This  is  shown  by  the  contours  in 
Fig.  16  where  W  at  the  corner  is  elevated  nearly 
five  orders  of  magnitude.  Refer  to  the  numerical 
values  in  Table  6. 

Localization  of  energy  density  at  the 
ceramic/rubber  interface  is  also  clearly  evi¬ 
denced  from  the  and  V  contours  in  Figs.  17 
and  18.  Fligh  values  of  S?  and  W  are  detected  at 
the  lower  side  of  the  ceramic  layer  where  it 
comes  into  contact  with  the  rubber.  This  effect  is 
particularly  pronounced  at  the  corner  of  elements 
22  and  25.  Concentration  of  at  the  top  of  the 
ceramic  layer  where  elements  23  and  24  meet  is 
also  noticeable  in  Fig.  17,  but  not  for  W  in  Fig. 


Table  6 

Volume  isoenergy  density  values  for  contours  in  Fig.  16 


Contour  'V  Contour  W 

No.  (Pa)  No.  (Pa) 


0.56X10  7  11 

0.76x10  4  12 


O.lSxHr-'  13 

4  0.23  14 

5  0.30  15 

o  0.38  16 

7  0.46  17 

8  0.53  18 

>1  0.61  19 

10  0.69  20 


0.76X10-1 

0.84 


0.99 

0.11  X10“2 
0.11 
0.12 
0.13 
0.14 
0.14 
0.15 


Corner  of 


Fig.  17.  Surface  energy  density  in  ceramic  layer. 


On 

1  <f  O 

[  Corner  of 

/  element  22  &  25 

1  2  2  / 

l _ _  1 _ s: 

If  1  r:  1 

Fig.  18.  Volume  energy  density  in  ceramic  layer. 
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Table  7 

Surface  isoenergy  density 

values  for  contours  in  Fig.  17 

Table  8 
Surface  iso 

energy  density  values  for  cont 

ours  in  Fig.  18 

Contour 

No. 

S’ 

(N/m) 

Contour 

No. 

S’ 

(N/m) 

Contour 

No. 

sr 

(Pa) 

Contour 

No. 

K- 

(Pa) 

1 

0.25X10-5 

11 

0.57X10“  2 

1 

0.53  X10“6 

11 

0.71X10 

2 

0.57X  10“3 

12 

0.63 

2 

0.71X10° 

12 

0.78 

3 

0.11X10"2 

13 

0.69 

3 

0.14X10 

13 

0.85 

4 

0.17 

14 

0.74 

4 

0.21 

14 

0.92 

5 

0.23 

15 

0.80 

5 

0.28 

15 

0.99 

6 

0.29 

16 

0.86 

6 

0.35 

16 

0.11X102 

7 

0.34 

17 

0.91 

7 

0.43 

17 

0.11 

8 

0.40 

18 

0.97 

8 

0.50 

18 

0.12 

9 

0.46 

19 

0.10  X 10“ 1 

9 

0.57 

19 

0.13 

10 

0.51 

20 

0.11 

10 

0.64 

20 

0.13 

21 

0.11 

21 

0.14 

18.  The  corresponding  contour  values  of  S?  and  Fig.  19  displays  the  S?  contours  which  tend  to 

W  can  be  found  in  Tables  7  and  8,  respectively.  congregate  at  the  corner  of  elements  21  to  26  at 
In  contrast  to  the  energy  density  distributions  the  top  and  elements  19  and  28  at  the  bottom, 

for  the  ceramic  layer  in  Figs.  17  and  18,  Figs.  19  Their  locations  with  reference  to  the  entire  lami- 

and  20  reveal  that  the  volume  energy  tends  to  nate  are  shown  in  Figs.  11  and  12.  The  highly 

dominate  in  the  rubber  rather  than  the  surface  packed  contours  of  W  near  the  top  and  bottom 

energy  density.  layer  in  Fig.  20  are  evidence  of  volume  energy 
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Table  9  Table  10 


Surface  isoenergy  density  values  for  contours  in  Fig.  19  Volume  isoenergy  density  values  for  contours  in  Fig.  20 


Contour 

No. 

(N/m) 

Contour 

No. 

(N/m) 

Contour 

No. 

V 

(Pa) 

Contour 

No. 

ar 

(Pa) 

1 

0.74x10  B 

11 

0.43X10  1 

1 

0.28x10”’ 

11 

0.88x10 

2 

0.43x10"- 

12 

0.47 

0.90X10° 

12 

0.97 

3 

0.86 

13 

0.52 

0.18X10 

13 

0.11  xio2 

4 

0.13x10  1 

14 

0.56 

4 

0.27 

14 

0.11 

5 

0.17 

15 

0.60 

5 

0.35 

15 

0.12 

6 

0.22 

16 

0.65 

h 

0.44 

16 

0.13 

7 

0.26 

17 

0.69 

7 

0.53 

17 

0.14 

8 

0.30 

18 

0.73 

8 

0.62 

18 

0.15 

9 

0.34 

19 

0.77 

9 

0.70 

19 

0.16 

10 

0.39 

20 

0.82 

10 

0.79 

20 

0.17 

21 

0.86 

21 

0.18 

trapped  at  the  top  interface  with  ceramic  and 
bottom  with  steel.  Tables  9  and  10  gives  the 
values  of  y  and  W. 

The  surface  and  volume  energy  density  con¬ 
tours  in  the  back-up  4340  steel  layer  are  shown  in 
Figs.  21  and  22.  Concentration  of  energy  density 
again  prevails  at  the  rubber /steel  interface  while 
the  energy  level  elsewhere  is  not  significant.  Nu¬ 
merical  values  of  the  contours  associated  with  y 


and  W  can  be  found  in  Tables  11  and  12,  respec¬ 
tively. 


5.  Conclusions 

Analysis  of  failure  behavior  in  composite  lami¬ 
nates  can  be  complex  if  detailed  information  is 
required  on  the  nonhomogeneous  nature  of  the 


, _ Corner  of 

'  element  21  &  26 


Fig.  20.  Volume  energy  density  in  rubber  layer. 
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Table  11 


Surface  isoenergy  density  values  for  contours  in  Fig.  21 

Contour 

S? 

Conte 

>ur  S? 

No. 

(N/m) 

No. 

(N/m) 

1 

0.51  XlO'3 

11 

0.71  X10'1 

2 

0.75  XlO'2 

12 

0.78 

3 

0.15  X 10"1 

13 

0.85 

4 

0.22 

14 

0.92 

5 

0.29 

15 

0.99 

6 

0.36 

16 

0.11X10° 

7 

0.43 

17 

0.11 

8 

0.50 

18 

0.12 

9 

0.57 

19 

0.13 

10 

0.64 

20 

0.13 

21 

0.14 

damage  process  that  changes  progressively  with 
the  load  history.  When  the  energy  is  transferred 
from  one  body  to  another  over  a  localized  region 
in  a  short  period  of  time  as  in  the  case  of  impact. 


Table  12 

Surface  isoenergy  density  values  for  contours  in  Fig.  22 
Contour  Contour  TtT 

No.  (Pa)  No.  (Pa) 

T  0.25X10“ 3  11  0.71X10 

2  0.71X10°  12  0.78 

3  0.14X10  13  0.85 

4  0.21  14  0.92 

5  0.28  15  0.99 

6  0.35  16  0.11  XlO2 

7  0.43  17  0.11 

8  0.50  18  0.12 

9  0.57  19  0.13 

10  0.64  20  0.13 

21  0.14 


additional  allowance  must  be  made  to  account  for 
the  wide  range  of  strain  and  strain  rates  that  the 
material  would  experience. 

Part  I  of  this  study  provided  the  methodology 


Fig.  22.  Volume  energy  density  in  4340  steel  layer. 
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for  deriving  the  dynamic  stress  and  strain  re¬ 
sponse  in  each  element  of  a  three-layered  com¬ 
posite  laminate  impacted  by  a  tungsten  rod  at 
high  speed.  The  concept  of  isoenergy  elements  is 
employed  such  that  the  base  material  properties 
obtained  from  uniaxial  tests  would  be  applied  to 
predict  the  behavior  of  elements  under  multi-axial 
stress  and  strain  states.  Time  scales  of  the  order 
of  10“'  p.s  were  found  to  coincide  with  the  pe¬ 
riod  of  initial  contact.  The  highly  nonhomogeous 
character  of  the  dynamic  disturbance  was  already 
evident.  The  surface  and  volume  energy  density 
near  the  impact  site  and  interfaces  were  elevated 
four  to  five  orders  of  magnitude.  Changes  in  the 
sign  of  the  stresses  are  indicative  of  the  effect  of 
wave  reflection.  Applicatioin  of  the  isoenergy 
density  theory  has  also  been  made  to  evaluate  the 
shock  response  of  undirectional  and  laminate 
composites  damaged  by  high  power  laser  [11]. 

Damage  evolution  of  the  ceramic/ rubber/ 
steel  laminate  is  presented  in  Part  II  [10].  As  the 
local  damage  alters  the  morphology  of  the  lami¬ 
nate  structure,  the  finite  element  grid  pattern 
needs  to  be  remeshed  for  each  time  step.  This 
process  is  continued  until  the  tungsten  rod  is 
completely  disintegrated  at  approximately  21.5  p,s 
after  impact.  Data  on  the  local  temperature  and 
disspation  energy  density  will  also  be  given  in 
Part  II  to  asertain  the  thermal  properties  of  the 
laminate. 
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